%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Fonction de calcul des "matrices de raideur" des efforts g�n�ralis�s
% �l�mentaires
%
% ne = nombre d'élements
% nse = nombre de sous-éléments par élement
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ Kf ] = Kf(poutre,sousPoutre,epsilon,a,b,nse,ne)

Kf1=zeros(4);
Kf2=zeros(4);

Kf1(1,1)=(9/4)*b*(1+2*a*epsilon);
Kf1(1,2)=(3/4)*b*(1-epsilon+4*a*epsilon);
Kf1(1,3)=-Kf1(1,1);
Kf1(1,4)=(3/4)*b*(1+epsilon+4*a*epsilon);

Kf1(2,1)=(3/4)*(6*epsilon*a*b+3*b-epsilon*b^2-2*a*(epsilon*a-1));
Kf1(2,2)=(1/4)*(12*epsilon*a*b+3*b*(1-epsilon)-epsilon*b^2-2*(1+epsilon*a)*(3*a-1));
Kf1(2,3)=-Kf1(2,1);
Kf1(2,4)=(1/4)*(12*epsilon*a*b+3*b*(1+epsilon)-epsilon*b^2-2*(1+epsilon*a)*(3*a+1));


for i=1:4
    
    Kf1(3,i)=Kf1(1,i);
    
end

Kf1(4,1)=(3/4)*(6*epsilon*a*b+3*b+epsilon*b^2-2*a*(epsilon*a+1));
Kf1(4,2)=(1/4)*(12*epsilon*a*b+3*b*(1-epsilon)+epsilon*b^2+2*(1+epsilon*a)*(3*a-1));
Kf1(4,3)=-Kf1(4,1);
Kf1(4,4)=(1/4)*(12*epsilon*a*b+3*b*(1+epsilon)+epsilon*b^2+2*(1+epsilon*a)*(3*a+1));


Kf1 = Kf1 * ((1/2)*(poutre.E*poutre.I)/((sousPoutre.L/nse)*(poutre.L/ne)^2));


for i=1:4
    for j=1:4
        if mod(i,2)==0
            Kf1(i,j)=(sousPoutre.L/nse)*Kf1(i,j);
        end
        if mod(j,2)==0
            Kf1(i,j)=(poutre.L/ne)*Kf1(i,j);
        end
    end
end







Kf2(1,1)=-(9/8)*a*epsilon;
Kf2(1,2)=-(3/8)*(3*a-1)*epsilon;
Kf2(1,3)=-Kf2(1,1);
Kf2(1,4)=-(3/8)*epsilon*(3*a-1);

Kf2(2,1)=(3/8)*(a-2*b)*epsilon;
Kf2(2,2)=(1/8)*(3*a-2*b-1)*epsilon;
Kf2(2,3)=-Kf2(2,1);
Kf2(2,4)=(1/8)*(3*a-2*b+1)*epsilon;


for i=1:4
    
    Kf2(3,i)=Kf2(1,i);
    
end

Kf2(4,1)=(3/8)*(a+2*b)*epsilon;
Kf2(4,2)=(1/8)*(3*a+2*b-1)*epsilon;
Kf2(4,3)=-Kf2(2,1);
Kf2(4,4)=(1/8)*(3*a+2*b+1)*epsilon;


Kf2 = Kf2 * ((poutre.E*poutre.I)/((poutre.L/ne)^3));


for i=1:4
    for j=1:4
        if mod(i,2)==0
            Kf2(i,j)=(sousPoutre.L/nse)*Kf2(i,j);
        end
        if mod(j,2)==0
            Kf2(i,j)=(poutre.L/ne)*Kf2(i,j);
        end
    end
end





Kf=Kf1+Kf2;


end